%Tobit regression model
function [theta]  = BetaTobit(rStock,rIndex, limit)
 
    [Hn Hu Hd] = ClassifyLimit(rIndex, rStock, limit);
    [beta, ~, r] = regress(rStock, rIndex);
    sigma = var(r);
    theta0 = [beta sigma];
    lb = [beta,sigma];
    ub = [beta*1.5,1];
    %Definitions Options
%     options = optimset('Algorithm','interior-point');
%     theta = fmincon(@BetaTruncatedLikelihood, theta0,[], [], [], [], lb, ub, [], options, Hn(:, 1), Hn(:, 2), theta0, limit);
    
    %Definitions Options
    options = optimset('Display', 'iter');
    options = optimset(options,'LargeScale','off');
    options = optimset(options,'GradObj','on','Jacobian','off');
    options = optimset('Algorithm','active-set');
    
%     if (isempty(Hu)==0 && isempty(Hd)==0)
%         theta = fmincon(@(x) -BetaLikelihoodTobit(x, Hn(:, 1), Hu(:, 1), Hd(:, 1), Hn(:, 2), limit), theta0, [], [], [], [], lb, ub, [], options);
%     elseif (isempty(Hu)==1 && isempty(Hd)==1)
%         theta = fmincon(@(x) -BetaLikelihoodTobit(x, Hn(:, 1), [], [], Hn(:, 2), limit), theta0, [], [], [], [], lb, ub, [], options);
%     elseif (isempty(Hu)==1)
%         theta = fmincon(@(x) -BetaLikelihoodTobit(x, Hn(:, 1), [], Hd(:, 1), Hn(:, 2), limit), theta0, [], [], [], [], lb, ub, [], options);
%     else
%         theta = fmincon(@(x) -BetaLikelihoodTobit(x, Hn(:, 1), Hu(:, 1), [], Hn(:, 2), limit), theta0, [], [], [], [], lb, ub, [], options);
%     end

    if (isempty(Hu)==0 && isempty(Hd)==0)
        theta = fminsearchbnd(@(x) -BetaLikelihoodTobit(x, Hn(:, 2), Hu(:, 1), Hd(:, 1), Hn(:, 1), limit), theta0, lb, ub);
    elseif (isempty(Hu)==1 && isempty(Hd)==1)
        theta = fminsearchbnd(@(x) -BetaLikelihoodTobit(x, Hn(:, 2), [], [], Hn(:, 1), limit), theta0, lb, ub);
    elseif (isempty(Hu)==1)
        theta = fminsearchbnd(@(x) -BetaLikelihoodTobit(x, Hn(:, 2), [], Hd(:, 1), Hn(:, 1), limit), theta0, lb, ub);
    else
        theta = fminsearchbnd(@(x) -BetaLikelihoodTobit(x, Hn(:, 2), Hu(:, 1), [], Hn(:, 1), limit), theta0, lb, ub);
    end
    
    theta = theta(1);
end